# coding:Utf-8

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

f = r'F:\Archive\newdata\aofeature.csv'
fb = r'F:\Archive\newdata\bofeature.csv'
a = r'F:\Archive\newdata\aopcafeature.csv'
b = r'F:\Archive\newdata\bopcafeature.csv'
tafile = r'F:\Resources\Train_A\Train_A\Train_A_001.csv'
adir = r'F:\Resources\Train_A\Train_A'
tdir = r'F:\Archive\TrainA'
data = pd.read_csv(f,header=None,usecols=list(range(0,60)))
datab = pd.read_csv(fb,header=None,usecols=list(range(0,60)))
pca = PCA(n_components=40)

j = pd.DataFrame(pca.fit_transform(data))
bj:np.ndarray = pca.transform(datab)
bj = pd.DataFrame(bj)
print(j)

"""
新特征是旧特征的线性组合，这些线性组合最大化样本方差，尽量使新的k个特征互不相关
说明矩阵在对应的特征向量上的方差越大，样本点越离散，越容易区分，信息量也就越多
pca 算法如下：
    对数据进行归一化处理（代码中并不是这么做的，而是直接减去均值）
    计算归一化后的数据集的协方差矩阵
    计算协方差矩阵的特征值和特征向量
    保留最重要的k个特征（通常k要小于n）。也能够自己制定。也能够选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k
    找出k个特征值相应的特征向量
    将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据
    
PCA 算法如下:
    组织数据形式，以便于模型使用；
    计算样本每个特征的平均值；
    每个样本数据减去该特征的平均值（归一化处理）；
    求协方差矩阵；
    找到协方差矩阵的特征值和特征向量；
    对特征值和特征向量重新排列（特征值从大到小排列）；
    对特征值求取累计贡献率；
    对累计贡献率按照某个特定比例选取特征向量集的子集合；
    对原始数据（第三步后）进行转换

"""
#计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征    
def meanX(dataX):
    return np.mean(dataX,axis=0)#axis=0表示依照列来求均值。假设输入list,则axis=1
"""
參数：
    - XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征    
    - k：表示取前k个特征值相应的特征向量
返回值：
    - finalData：參数一指的是返回的低维矩阵，相应于输入參数二
    - reconData：參数二相应的是移动坐标轴后的矩阵
"""
def pca__(XMat, k):
    average = meanX(XMat)
    m, n = np.shape(XMat)
    data_adjust = []
    avgs = np.tile(average, (m, 1))

    data_adjust = XMat - avgs

    covX = np.cov(data_adjust.T)   #计算协方差矩阵

    featValue, featVec=  np.linalg.eig(covX)  #求解协方差矩阵的特征值和特征向量 
    index = np.argsort(-featValue) #依照featValue进行从大到小排序
    finalData = []
    if k > n:
        print("k must lower than feature number")
        return
    else:
        #注意特征向量是列向量。而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值
        selectVec = np.matrix(featVec.T[index[:k]]) #所以这里须要进行转置 选取前k个重要特征向量
        finalData = data_adjust * selectVec.T  # 将原始数据投影到主成分上得到新的低维数据lowDDataMat
        reconData = (finalData * selectVec) + average  # 得到重构数据reconMat
    return finalData, reconData
finaldata,recondata = pca__(data,40)
print(finaldata)
# https://www.cnblogs.com/lzllovesyl/p/5235137.html
#coding=utf-8
from numpy import *

'''通过方差的百分比来计算将数据降到多少维是比较合适的，
函数传入的参数是特征值和百分比percentage，返回需要降到的维度数num'''
def eigValPct(eigVals,percentage):
    sortArray=sort(eigVals) #使用numpy中的sort()对特征值按照从小到大排序
    sortArray=sortArray[-1::-1] #特征值从大到小排序
    arraySum=sum(sortArray) #数据全部的方差arraySum
    tempSum=0
    num=0
    for i in sortArray:
        tempSum+=i
        num+=1
        if tempSum>=arraySum*percentage:
            return num

'''pca函数有两个参数，其中dataMat是已经转换成矩阵matrix形式的数据集，列表示特征；
其中的percentage表示取前多少个特征需要达到的方差占比，默认为0.9'''
def pca_(dataMat,percentage=0.9):
    meanVals=mean(dataMat,axis=0)  #对每一列求平均值，因为协方差的计算中需要减去均值
    meanRemoved=dataMat-meanVals
    covMat=cov(meanRemoved,rowvar=0)  #cov()计算方差

    U,sigma,VT = linalg.svd(covMat)



    eigVals,eigVects=linalg.eig(mat(covMat))  #利用numpy中寻找特征值和特征向量的模块linalg中的eig()方法
    k=eigValPct(eigVals,percentage) #要达到方差的百分比percentage，需要前k个向量
    eigValInd=argsort(eigVals)  #对特征值eigVals从小到大排序
    eigValInd=eigValInd[:-(k+1):-1] #从排好序的特征值，从后往前取k个，这样就实现了特征值的从大到小排列
    redEigVects=eigVects[:,eigValInd]   #返回排序后特征值对应的特征向量redEigVects（主成分）
    lowDDataMat=meanRemoved*redEigVects #将原始数据投影到主成分上得到新的低维数据lowDDataMat
    reconMat=(lowDDataMat*redEigVects.T)+meanVals   #得到重构数据reconMat
    return lowDDataMat,reconMat


#零均值化
def zeroMean(dataMat):      
    meanVal=np.mean(dataMat,axis=0)     #按列求均值，即求各个特征的均值
    newData=dataMat-meanVal
    return newData,meanVal
 
def percentage2n(eigVals,percentage):
    sortArray=np.sort(eigVals)   #升序
    sortArray=sortArray[-1::-1]  #逆转，即降序
    arraySum=sum(sortArray)
    tmpSum=0
    num=0
    for i in sortArray:
        tmpSum+=i
        num+=1
        if tmpSum>=arraySum*percentage:
            pass

def pca_P(dataMat,n):
    newData,meanVal=zeroMean(dataMat)
    covMat=np.cov(newData,rowvar=0)    #求协方差矩阵,return ndarray；若rowvar非0，一列代表一个样本，为0，一行代表一个样本
    eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特征值和特征向量,特征向量是按列放的，即一列代表一个特征向量
    eigValIndice=np.argsort(eigVals)            #对特征值从小到大排序
    n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n个特征值的下标
    n_eigVect=eigVects[:,n_eigValIndice]        #最大的n个特征值对应的特征向量
    lowDDataMat=newData*n_eigVect               #低维特征空间的数据
    reconMat=(lowDDataMat*n_eigVect.T)+meanVal  #重构数据

import numpy as np
from numpy.linalg import inv
from numpy import transpose as tr
# import ipdb


class PPCA(object):
    def __init__(self, q=2, sigma=1.0):
        self.q = q
        self.prior_sigma = sigma

    def fit(self, y, em=False):
        self.y = y
        self.p = y.shape[0]
        self.n = y.shape[1]
        if em:
            [self.w, self.mu, self.sigma] = self.__fit_em()
        else:
            [self.w, self.mu, self.sigma] = self.__fit_ml()

    def transform(self, y=None):
        if y is None:
            y = self.y
        [w, mu, sigma] = [self.w, self.mu, self.sigma]
        m = tr(w).dot(w) + sigma * np.eye(w.shape[1])
        m = inv(m)
        x = m.dot(tr(w)).dot(y - mu)
        return x

    def fit_transform(self, *args, **kwargs):
        self.fit(*args, **kwargs)
        return self.transform()

    def transform_infers(self, x=None, noise=False):
        if x is None:
            x = self.transform()
        [w, mu, sigma] = [self.w, self.mu, self.sigma]
        y = w.dot(x) + mu
        if noise:
            for i in range(y.shape[1]):
                e = np.random.normal(0, sigma, y.shape[0])
                y[:, i] += e
        return y

    def __ell(self, w, mu, sigma, norm=True):
        m = inv(tr(w).dot(w) + sigma * np.eye(w.shape[1]))
        mw = m.dot(tr(w))
        ll = 0.0
        for i in range(self.n):
            yi = self.y[:, i][:, np.newaxis] # np.newaxis==None
            yyi = yi - mu
            xi = mw.dot(yyi)
            xxi = sigma * m + xi.dot(tr(xi))
            ll += 0.5 * np.trace(xxi)
            if sigma > 1e-5:
                ll += (2 * sigma)**-1 * float(tr(yyi).dot(yyi))
                ll -= sigma**-1 * float(tr(xi).dot(tr(w)).dot(yyi))
                ll += (2 * sigma)**-1 * np.trace(tr(w).dot(w).dot(xxi))
        if sigma > 1e-5:
            ll += 0.5 * self.n * self.p * np.log(sigma)
        ll *= -1.0
        if norm:
            ll /= float(self.n)
        return ll

    def __fit_em(self, maxit=20):
        w = np.random.rand(self.p, self.q)
        mu = np.mean(self.y, 1)[:, np.newaxis]
        sigma = self.prior_sigma
        ll = self.__ell(w, mu, sigma)

        yy = self.y - mu
        s = self.n**-1 * yy.dot(tr(yy))
        for i in range(maxit):
            m = inv(tr(w).dot(w) + sigma * np.eye(self.q))
            t = inv(sigma * np.eye(self.q) + m.dot(tr(w)).dot(s).dot(w))
            w_new = s.dot(w).dot(t)
            sigma_new = self.p**-1 * np.trace(s - s.dot(w).dot(m).dot(tr(w_new)))
            ll_new = self.__ell(w_new, mu, sigma_new)
            print("{:3d}  {:.3f}".format(i + 1, ll_new))
            w = w_new
            sigma = sigma_new
            ll = ll_new
        return (w, mu, sigma)

    def __fit_ml(self):
        mu = np.mean(self.y, 1)[:, np.newaxis]
        [u, s, v] = np.linalg.svd(self.y - mu)
        if self.q > len(s):
            ss = np.zeros(self.q)
            ss[:len(s)] = s
        else:
            ss = s[:self.q]
        ss = np.sqrt(np.maximum(0, ss**2 - self.prior_sigma))
        w = u[:, :self.q].dot(np.diag(ss))
        if self.q < self.p:
            sigma = 1.0 / (self.p - self.q) * np.sum(s[self.q:]**2)
        else:
            sigma = 0.0
        return (w, mu, sigma)